差异基因没办法富集到通路你就放弃了吗

我在生信技能树分享了一个教程：不要怀疑，你的基因就是没办法富集到统计学显著的通路，强调了大家做生物信息学数据分析的同时，一定要加强统计学基础，比如把差异基因集（500个左右的基因）富集到KEGG数据库通路，本质上就是对每个通路单独做一次超几何分布检验罢了！

如果你的差异基因集（500个左右的基因）本身确实没有一些生物学功能的偏好性，非常有可能是在统计学显著阈值条件限定下，就是拿不到富集的kegg通路。这个时候，大家可以修改代码，比如：

kk.down <- enrichKEGG(gene         =  gene_down,
                        organism     = 'hsa', 
                        pvalueCutoff = 0.9,
                        qvalueCutoff =0.9)
# 需要自己差异分析筛选得到 gene_down 基因集 ，然后进行超几何分布检验

拿到通路后自己去调整阈值。

或者可以去MSigDB（Molecular Signatures Database）看看其它功能基因集是否可以富集到结果，MSigDB数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

• H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
• C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
• C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
• C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
• C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
• C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
• C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO  发表芯片数据
• C7: immunologic signatures: 免疫相关基因集合。

但是这个分析的前提是拿到差异基因集（500个左右的基因）

有些时候，大家挑选到的差异基因集数量过少或者过多，又不知道如何调整阈值，得到的结果会稀奇古怪，不可理喻。也就是说，这个分析严重依赖于差异基因集的挑选，因人而异，这样不好。

所以我们会推荐大家走GSEA分析，它并不会依赖于差异分析本身，我在生信技能树多次讲解GSEA分析：

在R里面的实现方式也是非常的简单：

# 需要自己对全部的基因进行一个排序，排序后的变量是geneList
library(clusterProfiler)
    kk_gse <- gseKEGG(geneList     = geneList,
                      organism     = 'hsa',
                      nPerm        = 1000,
                      minGSSize    = 10,
                      pvalueCutoff = 0.9,
                      verbose      = FALSE)
    tmp=kk_gse@result
    gsea_kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')

也可以对MSigDB数据库的全部基因集进行批量gsea分析，然后细致的去挑选结果。

# 选择gmt文件（MigDB中的全部基因集）
# 自己在 MigDB 官网下载 gmt 文件哦！
  d='~/biosoft/MSigDB/symbols/'
  gmts <- list.files(d,pattern = 'all')
  gmts
  #GSEA分析
  library(GSEABase) # BiocManager::install('GSEABase')
  ## 下面使用lapply循环读取每个gmt文件，并且进行GSEA分析
  ## 如果存在之前分析后保存的结果文件，就不需要重复进行GSEA分析。
  f='gsea_results.Rdata'
  if(!file.exists(f)){
    gsea_results <- lapply(gmts, function(gmtfile){
      # gmtfile=gmts[2]
      geneset <- read.gmt(file.path(d,gmtfile)) 
      print(paste0('Now process the ',gmtfile))
      egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
      head(egmt)
      # gseaplot(egmt, geneSetID = rownames(egmt[1,]))
      
      return(egmt)
    })
    # 上面的代码耗时，所以保存结果到本地文件
    save(gsea_results,file = f)
  }
  load(file = f)
  #提取gsea结果，熟悉这个对象
  gsea_results_list<- lapply(gsea_results, function(x){
    cat(paste(dim(x@result)),'\n')
    x@result
  })
  dev.off()
  gsea_results_df <- do.call(rbind, gsea_results_list)
   
  write.csv(gsea_results_df,file = 'gsea_results_df.csv')
  
  gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE') 
  gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6') 
  gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE') 
 

如果你差异基因没办法富集到通路，使用GSEA也许会有不一样的结果哦！

比如下面的这个通路，如果你使用上下调基因集去注释，通常呢你看不到它被富集到，因为top上下调基因都不是这个通路里面的，但是gsea分析就可以看到，它其实是显著下调的！